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ABSTRACT 

We present a method of reliably extracting the flux of individual sources from sky maps in the 
presence of noise and a source population in which number counts are a steeply falling function of 
flux. The method is an extension of a standard Bayesian procedure in the millimeter/submillimeter 
literature. As in the standard method, the prior applied to source flux measurements is derived from 
an estimate of the source counts as a function of flux, dN/dS. The key feature of the new method is 
that it enables reliable extraction of properties of individual sources, which previous methods in the 
literature do not. Wc first present the method for extracting individual source fluxes from data in 
a single observing band, then we extend the method to multiple bands, including prior information 
about the spectral behavior of the source population(s). The multi-band estimation technique is 
particularly relevant for classifying individual sources into populations according to their spectral 
behavior. We find that proper treatment of the correlated prior information between observing bands 
is key to avoiding significant biases in estimations of multi-band fluxes and spectral behavior, biases 
which lead to significant numbers of misclassified sources. We test the single- and multi-band versions 
of the method using simula ted observations wi th observing parameters similar to that of the South 
Pole Telescope data used in IVieira et al.l ()2009() . 

Subject headings: submillimeter — radio continuum: galaxies — methods: data analysis — methods: 
statistical 



1. INTRODUCTION 

The bias incurred when using noisy measurements to 
estimate source counts as a function of flu x has been 
a reco gnized issue since at l east the time of feddingtoiil 
()1913| ) . The seminal work of lScheuerl ()1957[ ) on the topic 
was crucial to reconciling apparently conflicting radio 
source count measurements and esta,blishi ng the model 
of an evolving universe (e.g., lLongaiilll966( ). 

Accounting for this bias becomes particularly impor- 
tant when the source population under investigation has 
counts that are a steeply falling function of flux, as is the 
case for sources in the 1-to-lOO mJy range at millime- 
ter/submillimeter (mm/submm) wavelengths. Sources 
in this flux and wavelength range have been of particular 
recent interest to astronomers, astrophysicists, and cos- 
mologists, primarily due to the recently disco vered popu- 
lation of high-rcdshift starburst galaxies (see iBlain et al.l 
()2002|) for a review), which has been shown to make up 
a substantial fraction of the cosmic infrared background 
(|Lagache et al.|[2005l ) and which provides strong tests for 
theories of galaxy and star formation. Both Bayesian 
and frcquentist methods have been develop ed for dealing 
with this bias in mrn/ subm m surveys (e.g.. lCoppin et aTl 



[200llMalonev et al.ll2005D . However, as will be discussed 
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in Sec. [51 the methods developed so far are not appropri- 
ate for estimating properties of individual sources. 

Dealing with this bias is further complicated when 
source properties ar e estimated f r om d ata in multiple 
wavelength bands. iMason et al.l (|2009[ ). for example, 
used a maximum-likelihood technique to make an unbi- 
ased estimate of the spectral index distribution of a pop- 
ulation of sources from two-band centimeter-wave data, 
but they did not attempt to estimate spectral proper- 
ties of individual sources. When multi-band data have 
been used to estimate properties of individual sources, 
the data in different bands have generally been treated as 
independent, despite the fact that the information used 
to correct the biased flux measurements is h ighly corre- 
lated between bands (e.g.. lGreve et al1l2008l) . As will be 
shown in Sec. [3l ignoring these correlations can result in 
severely misestimated source properties. 

In this work, wc propose a reliable, minimally biased 
estimator for the single-band brightness of individual 
sources. We then extend this formalism to the simul- 
taneous estimation of individual source properties over 
many bands, parameterized as either the source bright- 
ness in each band or as the brightness in one band 
and the spectral behavior between bands. The multi- 
band implementation of the method explicitly accounts 
for correlations in the information used to correct for 
flux bias in the individual bands. This work was moti- 
vated by the extragalactic sources detected in multi-band 
So uth Pole Tel e scope (SPT) data, which arc described 
in IVieira et al.l (|2009l hereafter V09); however, we be- 
lieve the method is directly applicable to other current 
or planned mm/submm experiments. We also believe 
the method should be appropriate for use in single- and 
multi-band surveys at other wavelengths, with the caveat 
that source variability could compromise the multi-band 
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version of the method if the observations in different 
bands are not simultaneous. 

2. SINGLE-BAND FLUX ESTIMATION 

The estimation of source brightness in mm/submm 
surveys is complicated by the fact that, at these wave- 
lengths, the number density of sources as a function of 
source flux is expected to be very steep. As a result, the 
measured flux of a source detected in a mm or submm 
survey will almost certainly suffer a positive bias, often 
referred to as "flux boosting." In this work, we define 
fl-irx boosting as the increased probability that a source 
we measure to have flux S is really a dimmer source plus 
a positive noise fluctuation over the probability that it is 
a brighter source plus a negative noise fluctuation. This 
asymmetric probability distribution means that naive es- 
timates of source fluxes will be biased high.^ In the litera- 
ture of mm/submm source surveys , one standar d method 
for dealing with this problem fe.g.. lCoppin et al . 2005) is 
to construct a posterior probability distribution for the 
intrinsi c flux of e ach d etect ion, as suggested by (a mong 
others) iJauiIcevI (I1968D and lHogg fc Turned (|1998D . Ac- 
cording to standard Bayesian reasoning, this posterior 
distribution is given by 



PiS,\S^) ^ P{S^\S.,)P{S,), 



(1) 



where Si is the intrinsic, true flux of the detected object, 
and Sm is the measured flux. P{Sm\Si) is the likelihood 
of measuring a flux Sm given a true flux Si, which in 
the simplest case is a Gaussian distribution centered on 
Si with width (T„, determined solely by the noise in the 
maps from which sources are being extracted. P{Si) is 
the prior probability of a source with intrinsic flux Si, 
which is proportional to the differential number counts 
vs. flux, dN/dS. 

However, the formulation in Eqn. [T] implicitly assumes 
that each detection corresponds to (at most) one real 
source. This is equivalent to assuming that the instru- 
ment used in the survey has infinite resolution. In real- 
ity, there will always be some possibility that a resolution 
element containing a source above the detection thresh- 
old will also contain one or more fainter sources which 
will contribute to the measured flux in that resolution 
element. As a consequence, even in a noiseless measure- 
ment, the probability of measuring a particular value of 
flux within a flnitc resolution element is not identical 
(in general) to the probability that a single source of 
that flux exists within the resolutio n element. In much 
of t he mm/submm literatur e (e.g., ICoppin et al.l l200"5l 
and lAustermann et al.l 120091) . P{Sp) (the probability of 
the total astrophysical flux Sp within a pixel) and P{S) 
(the probability of flnding a source of flux 5* within the 
solid angle of a pixel) are used interchangeably. In other 
words, what is being reported when this method is used 
is not the distribution of intrinsic source fluxes (which is 
a property solely of the source population) but rather a 

^ This phenomenon is closely rela ted to what is ref erred to in the 
literature as "Eddington bias" (e.g. . ITeerikorpill2004l ) : however, the 
consensus use of that term in the literature is to describe the bias 
introduced to estimation of source counts vs. brightness, not on the 
estimated brightness of individual sou rces. This usage is consistent 
with the original work of Eddington ( 19131 ) , and the distinction we 
dra w was also pointed ou t in (among others) rHoge fc Turner! II1998I ) 
and lCoppin et al.l ^OOEi). 



distribution of pixel fluxes, which is also dependent on 
the survey instrument. When calculating source counts 
using this method, it is possible to use simulations to ac- 
count for this instrumental depende nce or demonstrate 
that it has a negligible effect (e.g., lAustermann et ahl 
|2009() : however, the detailed shape of the posterior flux 
distribtutions for individual sources can still be affected, 
particularly at low signiflcance. 

One approach which avoids these particular compli- 
cations is to flnd the underlying (intrinsic) number 
counts model that is consisten t with the observed p ixel 
flux distributioii (as in, e.g.. iMalonev et all 120051 and 
iPatanchon et al.|[2009D or with the observ ed counts as a 
function of raw flux (as in "Reduction C" in lCoppin et al.l 
|2006|) . However, these methods are incapable of estimat- 
ing properties of individual sources. For this purpose, we 
extend the traditional Bayesian method to describe the 
properties of only one source in a given resolution ele- 
ment. Instead of defining the posterior with respect to 
the (intrinsic) total flux in a resolution element or pixel, 
we deflne the posterior with respect to the intrinsic flux 
of the single brightest source in the resolution element. 
The posterior probability distribution for this quantity 



P {SiyiaK\Sp^rn^ -^(iS'p^m | '5'max)-^(*5max) 5 
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where P(Sma.x\Sp^m) is the posterior probability that the 
true flux of the brightest source in a pixel is S'max given 
that we have measured the total flux in that pixel to be 
Sp.m] PiSpjn\Saiiix) is thc likclihood of measuring flux 
Sp.rn in a pixel given that the brightest source in that 
pixel has flux 5max; and -P( S'max) is the prior probability 
that the brightest source in that pixel has flux S'max- 

2.1. Prior Probability /'(S'max) 

The prior probability P(Smax) can be expressed as the 
probability that one source of flux Smax exists in that 
pixel times the probability that zero sources brighter 
than Smax exist in that pixel. As mentioned previously, 
the probability that within a pixel there exists a source 
of flux S is proportional to the differential number counts 
(per unit flux per unit solid angle) dN/dS. Under the 
assumption of purely Poisson statistics, the probability 
that zero sources above Smax exist in a pixel is 



P{N > Sn 



0) 
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(3) 



where fi is the mean number of sources above Sn 
pixels of size Afipi 

dN 



/J'(S'max) — Afij; 



dS 



dS. 



(4) 



In other words, P(S'max) should look like dN/dS with 
an exponential suppression at low S (where /i(S') will be 
large) : 



/"(S'max) OC 



dN 
dS 



exp 



s=s„ 



dN 
dS' 



dS' 



(5) 

We note that this choice of prior is natural in the con- 
text of characterizing individual source properties and 
avoids certain complications associated with the stan- 
dard choice of prior in the literature, which is the pixel 
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probability distribution that would be obtained from an- 
alyzing noiseless observations of a sk y with the assumed 
underlying source distri bution (e.g.. ICoppin erall[2005t 
lAustermann et all |2009[ ) . In particular, because point- 
source analyses of mm/sub- mm data almost inevitably 
involve spatially high-pass-filtering the data to remove 
large-scale noise and astronomical signals, the standard 
prior can take on negative flux values. This calls into 
question how such a prior can be describing the proba- 
bility of the intrinsic flux of an astrophysical object. 

2.2. Likelihood P{Sp^m\Si-niix) 

The total flux in a pixel given that the brightest source 
in that pixel has flux iSmax will be a sum of three contri- 
butions: the source at S'max, a- contribution from instru- 
mental or atmospheric noise in the survey, and a contri- 
bution from sources fainter than S'max-* We can think of 
these each having its own probability distribution, and 
the probability of the sum of their contributions to the 
flux in a pixel will be the convolution of the individual 
distributions (by the addition theorem for probabilities 
P{u ^ X + y) = J Px{x)Py{u — x)dx). We will assume 
that the contribution from instrument noise and atmo- 
sphere is well approximated by a Gaussian distribution: 



PiSp^m, noise-only) 



1 



(6) 



where cr„ is the width of the combined instrument and 
atmosphere noise distribution. The contribution from 
sources fainter than S'max is given by the probability of to- 
tal flux within a pixel, knowing that there are no sources 
in that pixel brighter than S'max- The probability for flux 
in a pixel given dN/ dS — the so-ca lled "deflexion p roba- 
bility" or P{D) — is worked out in IScheueil ()1957t ). and, 
under the assumption of pure Poisson statistics is 

P(5'p,™, noise-free) = FT{e^'''^^^"-'^°^^, (7) 

where FT{} denotes Fourier Transform, and r{uj) is the 
characteristic function of the probability of finding a 
source of flux S* in a pixel of solid angle AJl^: 



(8) 



To calculate the contribution of sources fainter than Smax 
under the constraint that there are no sources in a pixel 
greater than S'max, we apply the above result using a 
dN/dS distribution that is truncated at Smax. 

To summarize, the flux in a pixel given that we have ex- 
actly one source of flux Smax and zero sources above that 
flux will be be sum of the contribution from the source 
at Smax J the contribution from sources fainter than Smaxj 
and the contribution from noise. The probability distri- 
bution for this sum is the convolution of the individual 
distributions, so that 

ISmax) = (9) 

1 



.rn 1 5'niax) 



<5(Smax)*me[''^"'-''^°^l}* 



where "*" denotes convolution. 



® In reality, there will also be contributions from diffuse signals 
such as galactic dust or primary CMB anisotropy. However, these 
diffuse contributions will likely have been filtered out of maps used 
for source detection, so we neglect them here. 



2.3. Beam and filtering effects 

The exact formulations in Eqns. [S] and [5] only hold 
for a survey in which the only spatial filtering done in- 
volves binning into pixels. For most real instruments, 
the situation is complicated by the instrument beam (or 
point-spread function) and by the time-domain and map- 
domain filtering performed on the data. The likelihood 
in Eqn. IH] needs to be modified to refiect the difference 
in how sources fainter than Smax contribute to a real 
instrumental b eam cornpared to the top-hat pixel con- 
sidered earlier. IScheuerl ()1957l ) shows that Eqn. [8] can be 
modified for finite-resolution experiments by defining: 



r{uj) = FT 



dN 



dfl 



dShcam B ( 



(10) 



where B{9, (j)) is the angular response pattern of the in- 
strument, and 

S 

" B^M) ■ ^^^^ 
This formulation can be extended to account for any fil- 
tering in the data analysis by modifying the angular re- 
sponse function to include the filtering. 

The prior in Eqn. [5] and the (5(Smax) term in Eqn. [9] 
must also be modified for a finite-resolution experiment. 
We choose to define Smax for a finite-resolution exper- 
iment as the source that contributed the most fiux in 
a resolution element. That means that both equations 
must be replaced by integrals over the beam with Smax 
replaced by Smax times the beam response. Also, for 
each value of the integrand in Eqn. [101 the value at which 
dN/dS is truncated in the r{uj) calculation will be differ- 
ent. This quickly becomes computationally intractable, 
but here we are actually helped by the expected steep- 
ness of the source distribution. If we define Smax.p as 
the largest flux contribution in a resolution element, the 
probability that the source contributing that flux is a 
source of intrinsic flux very close to Smax.p that lies near 
the center of the beam is far larger than the probability 
that the contributing source is a much brighter source far 
off beam center. This means we can approximate the cor- 
rect versions of Eqs. [5] and [5] (that include integrals over 
the beam) with the original, inflnite-resolution versions, 
using a value for AJ7p that is roughly the area of the beam 
over which its response is near unity. In comparing to 
the simulated observations described in Sec. 12. 4[ we use 
1 arcmin^ , which is roughly the square of the full width 
at half maximum of the beam used in the simulations. 

2.4. Comparison with simulations 

We use simulated observations of mock skies includ- 
ing point-source populations drawn from model dN/dS 
distributions and Gaussian-distributed noise to test the 
formalism developed in the previous sections. Simulated 
observations were performed using three sets of observ- 
ing parameters, each roughly consistent with the 2.0 mm 
maps presented in V09, which we use as a concrete exam- 
ple. The three sets of observing parameters used are: 1) 
No spatial flltering beyond binning into 1-arcmin pix- 
els (and subtracting off the mean value of the map), 
noise consistent with the 2.0 mm SPT map shown in V09 
(~ 1.4 mJy rms); 2) Noise as in 1, but with spatial fllter- 
ing similar to the real SPT 2.0 mm maps in V09; 3) Fil- 
tering as in 1, but with the noise level halved. The source 
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count model that was used to cre a te the simulated maps 
is the sum of the iNegrello et all ([20?)1 850 /im counts 
for dust- dominated sour c es (sc aled to 2.0 mm as in V09) 
and the iDe Zotti et al.l ()2005[) counts for synchrotron- 
dominated sources. 

From these simulated observations, we extract the 
true, underlying posterior flux PDF, P{Sniax\Sp^rn), for 
three different values of Sp^m in each observing configura- 
tion. The true, underlying posterior flux PDF for a given 
value of Sp_rn was estimated by taking each source de- 
tected in the simulated maps with measured flux within 
SS of Sp^m, finding every source in the true, underly- 
ing source population that was associated with that de- 
tection, and recording the flux of the brightest associ- 
ated source as S'max for that detection. In the pixel-only 
cases, sources were considered associated with a detec- 
tion if they were in the same pixel as the detection; in 
the beam-and-filtering case, they were considered associ- 
ated if they were within 1 arcmin of the position of the 
detection. The true, underlying P(5'inax|'S'p^m) for that 
value of Sp_m is then simply the histogram of the S'max 
values assigned to all the detections with measured flux 
within SS of Sp^m- A total of 200 simulated 100 deg^ 
maps were used to construct the histogram. 

Fig. [T] shows the posterior flux PDFs extracted from 
the simulated observations and the calculated values of 
those posterior flux PDFs (using the formalism developed 
in the previous sections). The source model used to cal- 
culate the prior probability -P(5'max) is the same model 
used to generate the mock point-source skies. The values 
of Sp^rri for which posterior flux PDFs are shown corre- 
spond to detection significance values of 4.5(t, 5.5cr, and 
6.5a in each set of simulated observations. These detec- 
tion significance values correspond to raw flux values of 
6.3, 7.7, and 9.1 mJy in the full-noise simulations and 3.4, 
4.1, and 4.9 mJy in the halved-noise simulations. (The 
detection significance in the halved-noise simulations for 
a given raw flux value are not exactly twice those in the 
full- noise simulations because of the contribution of back- 
ground sources to the map rms.) 

For comparison with the output of the pixel-only simu- 
lated observations, the posterior flux PDF was calculated 
exactly using the equations in Sec. 12. H and Sec. 12.21 The 
top left and bottom left panels of Fig. [1] show that the 
calculated posterior PDF values for these cases are con- 
sistent with the simulation output to within the assumed 
Poisson errors. In the full-beam-and-flltcring case, r{uj) 
was calculated exactly, but the approximation described 
in the previous section was used for the prior and for the 
(^(Smax) term in Eqn.[9l The small but statistically mea- 
surable discrepancies between the simulated and calcu- 
lated PDFs in the fuU-beams-and-flltering case (middle 
left panel of Fig. [IJ are due to the imperfect nature of 
this approximation. 

2.5. Gaussian likelihood approximation 

The most computationally intensive step in estimating 
the posterior flux probability distribution P{S^a,y^\Sp^rn) 
is calculating the contribution from the other sources in 
the resolution element. For deep, single-dish mm/submm 
surveys, in which the angular resolution varies from ^ 10 
arcseconds to a couple of arcminutcs, this contribution 
will be dominated by the background of sources below 
the confusion limit (the regime in which there are many 



sources per resolution element). By the central limit the- 
orem, as the number of sources per resolution element 
becomes large, the distribution of total flux from these 
sources will approach a Gaussian. If we use this fact to 
approximate the background source contribution to the 
posterior as another Gaussian-distributed noise source, 
the calculation of the likelihood in Eqn. [5] becomes triv- 
ial: 

P(S \S ) - ^ p-(Sp.„-5„a.-s;)V2at"ot 

-'^ \^p,m\ ^ max) — > 7^ — ^ 1 

- . . . . 

where Sp is the mean contribution from sources in the 

map, and atot is the quadrature sum of the instrumental 
and atmospheric noise and the fluctuations in the source 
background: 

There will always be a high-flux tail to the P{D) distri- 
bution, because some resolution elements in the survey 
will have more than one source brighter than the con- 
fusion limit. However, for steep source populations, the 
mean number of sources contributing to the total flux 
in a pixel will be large, and the tails will be less im- 
portant. Furthermore, because the background source 
contribution adds in quadrature to the instrument and 
atmospheric noise contribution (which is likely to be very 
well approximated by a Gaussian), non-Gaussianity in 
the source contribution will manifest itself in the total 
P{D) distribution only if it is the dominant source of 
noise in a survey. For a survey like the SPT, this is 
not the case, and the Gaussian likelihood approximation 
holds very well, as shown in the top and middle right 
panels of Fig. [TJ However, a survey that was a factor of 
two deeper at the same wavelength and resolution would 
incur much more significant errors by adopting the Gaus- 
sian likelihood approximation, as shown in the bottom 
right panel of Fig. [TJ 

It is worthwhile to note that adopting the Gaussian 
likelihood approximation removes one of the differences 
between the fiux estimation method presented here and 
the standard method in the literature. In the Gaus- 
sian likelihood approximation, the implicit assumption 
is made that the only contributions to the measured flux 
in a pixel are instrumental and atmospheric noise, a sin- 
gle bright source, and a source background that acts as 
another symmetric noise source. In this case, the only 
difference between the likelihood of measuring some flux 
in a pixel given the flux of the brightest source in that 
pixel and the likelihood of measuring some flux in a pixel 
given the intrinsic total flux in that pixel is the mean 
contribution from the source background. Since all of 
the mm/submm experiments under consideration here 
are differential (i.e., not sensitive to the mean brightness 
on the sky), the two likelihoods are effectively identical. 
However, a significant difference remains between the pri- 
ors used in the two methods, as discussed in Sec. 12.11 

2.6. Estimating source counts with a source-count prior 

If the individual source fluxes estimated with this 
method (or with the standard method in the literature) 
are used in estimating source counts (as they are in V09 
and most submm analyses), one might object that a prior 
probability has been assumed for the very quantity that 
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Fig. 1. — Left Panels: True, underlying posterior flux PDF P(Smax|S'p.7n) (extracted from simulated observations and shown by symbols with 
error bars) and calculated values for that PDF (using the procedure outlined in Scc.[2]and shown by lines). In the top two rows, the raw flux values 
for which the posterior PDFs are calculated are 6.3 mjy {diamond symbols and dashed line), 7.7 mjy {triangle symbols and dot-dashed line), and 
9.1 mJy {square symbols and triple dot-dashed line). In the bottom row, those raw flux values are 3.4, 4.1, and 4.9 mJy. In both cases, these raw 
flux values correspond to detection significance levels of 4.5cr, 5.5tT, and 6.5cr. Vertical error bars on the extracted PDFs are from Poisson statistics; 
horizontal error bars show SS — 0.2 mJy, the width of the flux bins used to construct the -P(Smaxl'S'p.m) histogram from the simulated observations 
(see text for details). Right Panels: Gaussian likelihood approximation to the posterior PDF calculation. Solid lines are from the full calculation 
and are identical to the lines in the left panel; dashed lines are calculated using the Gaussian approximation to the likelihood P(S'p. m I'S'max) shown 
in Eqn. [12] Top R,ow: binning into 1-arcmin pixels only, 1.4 mJy noise rms; Middle R,ow: full beam and filtering, 1.4 mJy noise rms; Bottom 
Row: binning into 1-arcmin pixels only, 0.7 mJy noise rms; 
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one is attempting to measure. It is important to remem- 
ber that the source-count prior is apphed probabihsti- 
cally to determine each source's flux, so the estimated 
source counts will only resemble the prior if the measure- 
ments have no constraining power. From the posterior 
flux PDFs shown in Fig. [U it is clear that for the combi- 
nation of data and prior used in these simulated observa- 
tions, the data are providing the bulk of the information 
in the flux measurement down to the lowest detection 
significance shown (4.5cr). For a real survey, plots of 
individual posterior flux PDFs such as these are useful 
for determining where the reported source counts contain 
significant new information and where they simply repro- 
duce the prior. Another useful check on the constraining 
power of the data is to vary the source count prior and 
confirm that the estimated cou nts do not change signifi- 
cantly (as in lScott et alll2008l and V09). 

3. MULTI-BAND SOURCE FLUX ESTIMATION 

The posterior fiux estimation method described in 
Sec.[2]implicitly assumes that both the data and the prior 
knowledge of the source distribution are restricted to a 
single observing band. The situation will inevitably be 
more complicated when both data and priors are avail- 
able in multiple bands. In full generality, the task at 
hand is now to simultaneously estimate the posterior 
probability of the intrinsic source flux in multiple bands 
given the measured flux in those bands and any prior 
information. Using the formalism from Sec. [5J we can 
write the two-band case as 

) cx (14) 

In the mathematically simplest case, both the condi- 
tional and prior probabilities on the right-hand side of 
Eqn.[T4]are uncorrelatcd between bands, so that the pos- 
terior probability distribution for flux in the two bands is 
simply the product of the two single-band distributions: 

-P(*5'inax,l 1 "5'max.2 I'S'p^m.l 7 '^Pi'n.z) OC (15) 
,m,l I «5'max,l )P(^max.l)P(S'p ,m,2 I •5'max,2 )P{S max, 2 ) • 

However, the assumption of uncorrelatcd prior informa- 
tion in the two bands will very rarely be valid. The prior 
in each band will be derived from previous estimates of 
source counts at wavelengths and flux levels as near as 
possible to that band, using assumptions about spectral 
behavior to extrapolate to the bands of interest. For 
bands that are reasonably close to each other in wave- 
length, the existing data that go into estimating the pri- 
ors in the two bands will almost certainly have some over- 
lap. This is not simply a matter of imperfect priors; in 
the limit of perfect prior knowledge of the source counts 
in both bands, the priors in the two bands would only be 
independent if the source populations measured in the 
two bands had zero overlap. In the example of the SPT 
data at 1.4 and 2.0 mm analyzed in V09, the prior prob- 
ability P(5'max) in both bands is dominated in the few- 
mJy flux region by the dusty starburst population, with 
some contribution from synchrotron-dominated AGN. It 
is clearly a poor approximation to assume that the priors 
in these two bands are uncorrelatcd. 

It is possible to construct the full two-dimensional prior 
probability -P(5'max.i, »S'max.2); however, it will often be 



more convenient to change variables and cast this prior 
probability in terms of the flux in one band and the ex- 
pected spectral behavior of the sources contributing to 
the prior. We define the spectral index a through the 
relation 

5(A2)-5(Ai) " (16) 

and write the two-dimensional prior as /'(S'max,!, ot). The 
only requirement for this prior to factor into independent 
priors P(5max,i) and P{a) is that the spectral index dis- 
tribution between the two bands not depend on flux. Of 
course, in reality P{a) will depend somewhat on flux — 
for example, in the SPT case, the brightest sources are 
mostly AGN, but near the confusion limit the source pop- 
ulation is dominated by dusty galaxies. In cases such as 
the analysis of SPT data in V09, in which estimating the 
spectral index distribution of the sources is one of the key 
goals of the analysis, the prior used will be sufficiently 
weak (V09 use a flat prior from —3 < a < 5) that the 
slight dependence on flux of the real P{a) is of no con- 
sequence. Alternatively, we can accept the added com- 
plexity of using the full two-dimensional prior on Smax 
and a — assuming sufficient data exist to meaningfully 
construct such a distribution. 

No matter how we choose to construct the spectral 
index prior, we now have the two-dimensional posterior 
PDF of flux in one band and a: 

P{S 

max, 1 ! Q;| 5*^^772^1 , 5*^.771^2 ) (1'^) 
P{Sp^m,l, Sp^m,2\Smax,l, Q^)-P(5'max,l , Q^)- 

If we choose, we can then transform this into a two- 
dimensional posterior probability distribution for the flux 
in both bands: 



-P(»S'max,lj 5niax,2 |>S'p,m,l I 5'p, 771,2) 



(18) 



-P('S'max,lj Q^jiSp, 771,1, Sp,nn,2) 



da 



dSniax,2 

where <ia/dS'max,2 is derived from Fan. [T51 

3.1. Choice of detection band 

We note that as the prior on a is relaxed to a flat 
prior between — oo and oo (independent of source flux), 
the posterior PDF for S'max.i in Eqn. [T71 reduces to the 
single-band posterior of Eqn. [21 Meanwhile, the pos- 
terior PDF for S'max.2 bccomcs equal to the likelihood 
P(5'p, 771, 2 1 S'max, 2) — i.e., there is effectively no prior on 
<S'max,2- This poiuts out an apparent asymmetry between 
bands in our approach, namely that the two-band pos- 
terior flux PDF will depend on which band you choose 
to use prior source count information from — call this 
the "detection band" — and which band you apply prior 
information to only through the combination of the prior 
on a and source-count prior on the first band. 

The real issue is that dN/dS in each band and a dis- 
tribution in a are not three independent pieces of infor- 
mation; rather, the choice of dN/dS in one band and a 
distribution in a uniquely specifies dN/dS in the other 
band. If the assumptions regarding dN/dS in both bands 
and the distribution of a arc internally consistent, then 
the two-band posterior fiux PDF will be identical regard- 
less of which band is chosen as the detection band. How- 
ever, if one of the goals of the analysis is to measure the 
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distribution of spectral indices, such a prior on a may 
be too restrictive. In V09, the choice was made to adopt 
a non-restrictive prior (—3 < a < 5) and analyze the 
data using each band in turn as the detection band. The 
differences in the final derived source counts between the 
two analyses were small compared to the statistical un- 
certainties. 

3.2. Variable sources 

The members of at least one of the populations ex- 
pected to contribute to mm/submm source counts — 
namely AGN — are expected to have highly time- 
variable flux. This complicates any estimate of source 
properties from data taken at different epochs. For in- 
struments that observe simultaneously in multiple wave- 
length bands — as the SPT does — this is not an is- 
sue, and we have ignored any effects of variability on the 
method presented here. This would be an issue, however, 
in extending this method to radio survey data. 

3.3. Gaussian likelihood approximation in the 
multi-band case 

As in the single-band case, the calculation of the like- 
lihood P{Sp^rn.i, Sp^m.2\^nia.x,i,a) IS greatly simplified by 
the assumption that the instrumental and atmospheric 
noise and the contribution from sources below Smax are 
Gaussian-distributed. While the full, non-Gaussian two- 
dimensional likelihoo d is in pr i nciple calculable — by ex- 
tending the result of iScheueil ()1957[ ) into a multivariate 
Fourier operation — that derivation is beyond the scope 
of this work. It is also possible to estimate the full, non- 
Gaussian two-dimensional distribution by simulated ob- 
servations. For now, we will adopt the Gaussian likeli- 
hood approximation and use simulated observations to 
check its validity, as we did for the single-band results. 

Under the Gaussian likelihood approximation, the two- 
band likelihood (analogous to the single-band likelihood 
in Eqn. I12p is given by 



P{Sp.m,l: *5'p,m,2 I'S'max.l I Oi) 



exp (-i r^C-ir) 



(19) 



27rVdctC 

where C is the noise covariance between the bands 
(including contributions from instrument noise, atmo- 
sphere, and sources fainter than S'max), and r is the resid- 
ual vector 



1^ 1 *5'niax,l *^p,l? 
5'p,m,2 — 'S'max,2(Q!) ^ 'S'p,2| 
■ 771^1 »5'niax,l 



(20) 



5. 



p,m,2 



Smax,l 1 ) "^^,2 



3.4. Comparison to simulations 

Analogous to Fig. [T] for the single-band case. 
Fig. [5] shows comparisons between calculated ver- 
sions of the two-dimensional posterior distribution 
^'(S'max,i,'S'max,2|5'p,m,i,5'p,,„,2) and the true values of 
those distributions (extracted from simulated observa- 
tions). The simulated observations were performed 
by populating fake skies at two observing wavelengths 



20 



15 



s 10 




8 10 



„[mJy] 



20 



^ 15 hi'. 

E 



S 10 




10 



„ [mJy] 



Fig. 2. — Comparison of the true, underlying posterior two-band 
flux PDF P(5'max,l , 5'max,2 |'S'p,TTi,l , 'S'p,TTi,2) (extracted from simulated 
observations) and calculated values for that PDF for two possible pairs 
of measured 1.4 mm and 2.0 mm flux. Top Panel: Sp^m — 14.2 mjy 
(4.2cr) at 1.4 mm, Sp,™ = 7.0 m,Jy (4.9a) at 2.0 mm. 2.0 mm used 
as "detection band" (see See. 13. ll ■ Bottom Panel: Sp^m — 17.2 mJy 
(4.9cr) at 1.4 mm, Sp,m = 6.0 mJy (4.1a') at 2.0 mm. 1.4 mm used 
as "detection band" . Symbols and contours are as follows: Black X: 
Value of measured flux Sp^m in the two bands. Black (solid) con- 
tours: true, underlying posterior two-band flux PDF, extracted from 
the simulated observations as described in Sec. 13.41 (These contours 
have been smoothed slightly.) Blue (dashed) contours: posterior PDF 
calculated using the multi-band source flux estimation procedure de- 
scribed in Sec.[3]with a flat prior on spectral index a with cutoff values 
of —3 and 5. Red (dotted) contours: similar to blue (dashed) contours, 
but with a prior on a equal to the spectral index distribution used in the 
input to the simulated observations (a — 2.7 it 0.3). Green (dash-dot) 
contours: posterior PDF obtained using the single-band procedure de- 
scribed in Sec.[2]independently in both bands, ignoring the correlations 
betwee n the priors in the tw o bands. All posterior PDF calculations 
use the lNegrello et al.l (I2007f ) source counts model to construct the flux 
prior. All contours are drawn at 0.5ct intervals. 

(1.4 mm and 2.0 mm) with a single population of 
sources with an underlying Gaussian distribution of spec- 
tral indices a = 2.7 ± 0.3, similar to the spectral in- 
de x distributio n for h igh-redshift, dusty galaxies derived 
in iKnox et al.l ()2004| ). The number counts as a func- 
tion of flux for the source population come from the 
iNegrello et al.l (|2007[ ) dusty galaxy counts at 850/.tm; this 
model was also used to construct the source flux prior. 
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Noise was added to the fake skies in each band at a 
level similar to that in the corresponding SPT bands in 
V09. To avoid confusing effects of slight misestimations 
of priors due to beam and filtering with fundamental is- 
sues in the two-band implementation, this set of sim- 
ulated observations involved no spatial filtering beyond 
binning into 1-arcmin pixels. Posterior two-band flux 
PDFs were extracted from the simulated observations as 
in the single-band case, namely by finding the brightest 
source associated with each detection in the simulated 
maps and constructing a two-dimensional histogram of 
{'5'max.i, '5'max,2} for every pair of measured flux values. 

True and calculated values of the posterior PDF 
are shown in Fig. [2] for two values of measured two- 
band flux. These raw flux values at 1.4 mm and 
2.0 mm — {14.2,7.0} mJy in the top panel of Fig. [5] 
and {17.2,6.0} mJy in the bottom panel of Fig. [5] — 
correspond to detection significances of {4.2,4.9} and 
{4.9,4.1}. These values were chosen to illustrate the 
importance of using the full two-band information as 
opposed to calculating each band's posterior flux PDF 
individually and ignoring correlations in the two bands' 
prior information. 

Three different calculated values of the posterior are 
shown in each panel. Two versions of the calculated 
posterior use the two-band implementation described in 
Sec. m one using a fiat prior on a (—3 < a < 5) and 
one using a prior on a that is equal to the true un- 
derlying a distribution (a Gaussian with a = 2.7 and 
(Tq = 0.3). In each of these cases, the band in which the 
source was detected more significantly was used as the 
detection band. The last version of the two-band pos- 
terior flux PDF shown in Fig. [5] is the product of the 
individual single-band posterior PDFs, each calculated 
using the procedure outlined in Sec. [5] and assuming no 
correlations between the prior information in each band. 

The large amount of information in Fig.[2]can be boiled 
down to three main points: 

1. If one had perfect prior information on dN /dS and 
a, one could calculate the posterior two-band flux 
PDF for every source perfectly. 

2. With far less restrictive priors on the a distribu- 
tion, one can make an estimate of the posterior two- 
band flux PDF for every source that has no strong 
bias but has somewhat less constraining power. 

3. Calculating the posterior two-band flux PDF by 
de-boosting each source individually and assuming 
no correlation between the priors in each band can 
result in highly biased posterior distributions. The 
posterior flux estimate in the band in which the 
source is detected strongly is reasonable, but the 
fiux estimate in the other band is de-boosted to 
the confusion limit. The spectral index inferred 
from this calculation is actually a worse estimate 
of the true index than using the measured fiux in 
each band uncorrected for boosting (as shown by 
the black crosses in Fig. [5]). This highlights the 
importance of accounting for the correlations in 
the prior information used to de-boost multi-band 
fiuxes. particularly when dealing with a source pop- 
ulation for which the number counts are a steeply 
falling function of fiux. 



3.5. Classification of sources through their spectral 
behavior 



q 



0.4 




i ' ' ' 1 1 1 1 

1 

\ 

\ G] : 


0.3 




\ ^ 


0.2 




\ ^ 

\ ^ ■ 

\ ' 
\ ^ 


0.1 






0.0 




~ 


4 


5 6 7 8 E 






2 mm detection significance 


0.4 




1 


0.3 




p : 
/ \ : 


0.2 


/ 

/ 


/ ^ 1 

1 ^ \ 
g : 
\ ; 

\ ' 


0.1 


/ 

/ 

; 


\ ~ 


0.0 


: «^ 


< < < < 1 1 1 1 



4 5 6 7 8 9 

2 mm detection significance 



Fig. 3. — Fraction of sources in simulated observations which would 
be misclassified by a posterior spectral index PDF criterion which di- 
vides sources into synchrotron- or dust-dominated populations based 
on the value of P(a > 1.5). See Sec. 13.5) for details of the simulated 
observation. Diamonds and solid lines show the misclassification frac- 
tion obtained using the two-band posterior flux estimation described 
in Sec. [3] boxes and dashed lines show the misclassification fraction 
obtained using the single-band procedure described in Sec. f2]indcpen- 
dently in both bands and ignoring the correlations in the two bands' 
priors. Top Panel: fraction of all sources misclassified as a function 
of 2.0 mm detection significance. Bottom Panel: fraction of sources 
misclassified which had P(a > 1.5) > 0.9 or P(a > 1.5) < 0.1, also 
as a function of 2.0 mm detection significance. The feature at ~ 4.7o' 
in the dashed curve in the bottom panel is an artifact of the specific 
noise levels chosen for the two bands (for details, see Sec. I3.5T . 

One key use of spectral information is to separate 
sources into different populations. For example, V09 use 
the posterior probability distribution of a for every de- 
tected source to classify it as either synchrotron- or dust- 
dominated. This allows the source counts in each popu- 
lation to be compared to models of that population and 
to counts at other wavelengths where that population is 
dominant. Any bias in the posterior two-band flux PDF 
of a source could cause sources to be misclassifled, lead- 
ing to a bias in the estimation of source counts for both 
populations. 
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We have used simulated observations very similar to 
those used in Sec. 13.41 to investigate how often sources 
are misclassified using the two-band posterior flux PDF 
estimated with and without accounting for correlations 
in the prior information between bands. The only dif- 
ference in this set of simulated observations and those of 
Sec. 13.41 is that there are two populations used to create 
the fake skies that are observed by our instrument. The 
two populations are the dust-dominated population used 
in Sec. l3.4l Dlus a sy nchrotron-dom i nated population with 
source counts as in lDe Zotti et al.l (|2005[ ) and a Gaussian 
spectral index distribution with a = — 0.7± 0.5, roughly 
consistent with the spectral behavior of the brightest 
synchrotron-dominated sources in V09. 

For each detection in the simulated maps, the brightest 
source associated with that detection was identified, and 
the posterior spectral index PDF for the source was cal- 
culated by initially calculating the two-dimensional pos- 
terior PDF for flux in one band and a and then marginal- 
izing over the flux variable to create a one-dimensional 
spectral index PDF. For this exercise, we only use the 
broad, flat a prior from Sec. 13.41 (—3 < a < 5). The 
posterior a. PDF was then compared to the true spectral 
index of the brightest source associated with that detec- 
tion. We classified detections in the simulated maps as 
synchrotron-dominated if P{a > 1.5) < 0.5 and dust- 
dominated if P{a > 1.5) > 0.5; similarly, we classi- 
fied the brightest source associated with the detection 
as synchrotron- or dust-dominated according to whether 
its true spectral index was greater than or less than 1.5.^ 
The posterior a PDF for each detection was also esti- 
mated by calculating the posterior flux PDF in each band 
independently using the procedure outlined in Sec. [2] — 
ignoring any correlations between the prior information 
in the two bands — and combining the two flux PDFs 
to create a PDF for a. This a PDF was used to classify 
sources similarly to the two-band a posterior. 

The fraction of sources misclassified (labeled as 
synchrotron-dominated using the posterior a PDF 
when the brightest associated source was in fact dust- 
dominated, or vice- versa) in each case is shown as a func- 
tion of single-band detection significance in Fig. [3l The 
two-band implementation — even with the weak prior on 
a — shows a clear improvement in misclassification frac- 
tion over combining independent single-band PDFs at all 
significance levels up to 7 a. Particularly striking is the 
difference in "high-confidence" misclassifications — in- 
stances in which the classification based on the posterior 
a PDF was at the 90% confidence level or greater, but 
was wrong. As shown in the bottom panel of Fig.[3l the a 
estimation based on independent single-band PDFs has 
up to a 25% rate of high-confidence misclassifications, 
but rate for the two-band implementation is effectively 
zero. 

The turnover at ~ 4.7(t in the rate of high-confidence 
misclassifications using the single-band information and 
ignoring correlations (Fig. [3l bottom panel, dashed 
curve) is due to the relative noise level in the two bands. 
Because the noise at 1.4 mm is more than twice that at 

^ The threshold value of a = 1.5 was chosen to lie roughly at 
the minimum of the two-population spectral index histogram (for 
sources above the detection threshold) The results in this section 
are insensitive to moving this threshold value by ±0.5. 



2.0 mm, there are many sources that are intrinsically 
dust-dominated which are nevertheless detected more 
significantly at 2.0 mm. If such a source is detected above 
4.5cr at 2.0 mm but below 4.5cr at 1.4 mm, the posterior 
fiux PDF for that source will be centered near the raw, 
detected fiux at 2.0 mm but de-boosted to the confusion 
limit at 1.4 mm. This results in a robust "measurement" 
of a negative spectral index for this source and, hence, 
a high-confidence misclassification. If the source is de- 
tected below 4.5(7 in both bands, there is effectively no 
constraint on a using this procedure, so the source may 
be misclassified, but not at high confidence. 

3.6. More than two bands 

The two-band formalism laid out in Sec. |3] is suffi- 
cient for the V09 analysis of two-band SPT data. How- 
ever, data in three or more bands at mm/submm wave- 
lengths and mjy fiux levels have been or are currently be- 
ing collected by the Balloon-bor ne Large- Aperture Sub- 
miUimeter Telescope^° (BLAST. IPevlin et al.l[2009h . the 
SPT, and the Ata cama Cosmology Telescope^^ (ACT, 
iFowler et all 120071 ). Further more, the recent launch of 
Herschel"'^^ and Planck-'^'^ means that we will soon have 
simultaneous measurements of mm/submm sources in as 
many as seven bands (depending on where you choose to 
define the limits of the mm/submm spectral region). 

Fortunately, the two-band formalism laid out in Sec. [3] 
is easily extended to more than two bands, although the 
calculation necessarily becomes more complex. For the 
case in which the Gaussian approximation holds and each 
source's spectral behavior can be described by a single 
power-law index across all bands, then the multi-band 
calculation is a trivial extension of Eqns.[l7]-[20l A first 
step in relaxing the assumption of a single spectral index 
for each source would be allowing a break in the spec- 
trum such that each source would have a single spec- 
tral index between any two bands. The spectral index 
prior would then be a function of A'bands — 1 variables 
a = {ai2, Q!i3, aiAf}. These variables would neces- 
sarily be highly correlated, so we would require the full 
{N — l)-dimensional prior. A full treatment of the A^- 
band version of the method, including tests on simulated 
observations, will be the subject of future work. 

4. CONCLUSIONS 

We have constructed a method for reliable, minimally 
biased estimation of single-band and multi-band prop- 
erties of individual sources from noisy data. We find 
that proper treatment of correlated prior information in 
the multi-band version of the method is crucial to avoid 
significant biases in estimates of multi-band fluxes and 
spectral behavior. The single- and multi-band imple- 
mentations of the method have been verified through 
simulated observations of mm data, and the two-band 
implementation has been used to est imate source fluxe s 
and spectral behavior in SPT data (jVieira et al.ll2009t ). 
This method, or an extension thereof to more than two 
bands, is directly applicable to source analyses for most 
current and upcoming mm/submm experiments, includ- 

http: //www.blastexperlment . In f o| 
http : //www . physics . pr inceton . edu/act/ | 
http: //www. esa. int/science/herschel 
13 http://www7esa.int/SPECIALS/Planck/ 
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ing BLAST, ACT, Planck, and Hcrschcl and should also 
be applicable to data taken at other wavelengths. 
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